Quantum Chaos in the Bose-Hubbard model 
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We present a numerical study of the spectral properties of the ID Bose-Hubbard model. Unlike 
the ID Hubbard model for fermions, this system is found to be non-integrable, and exhibits Wigner- 
Dyson spectral statistics under suitable conditions. 

PACS numbers: PACS: 03.65.-w, 32.80.Pj, 03.75.Nt, 71.35.Lk 



During the last two decades much attention has been 
paid to the spectral properties of quantum systems which 
exhibit chaotic dynamics in the classical limit. The in- 
trinsic complexity of these systems (due to underlying 
chaos) leads to an irregular energy spectrum with univer- 
sal statistical properties - a phenomenon often referred to 
as 'quantum chaos' P|. One of the recent developments 
in the field is related to 'many-body quantum chaos' 0- 
Here the systems of interest consist of many interacting 
(identical) particles and usually have no classical coun- 
terpart. Notwithstanding, numerical studies of some sim- 
ple models (like, for example, a ID chain of interacting 
spins (3|) have shown that the spectrum of a deterministic 
many-body system can also exhibit universal statis- 
tics. In the present paper we numerically study the spec- 
trum of the ID Bose-Hubbard model 5]. Our particular 
interest in this system is due to its possible laboratory re- 
alization with cold atoms loaded into optical lattices 
It is shown below that, unlike the ID Hubbard model for 
fermions, the ID Bose-Hubbard model is generally non- 
integrable. Moreover, in a certain parameter regime its 
spectrum obeys the Wigner-Dyson statistics for Gaus- 
sian Orthogonal Ensembles and, hence, this system can 
be considered as a quantum chaotic system. 

The Hamiltonian of the ID Bose-Hubbard (BH) model 
reads 
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Having in mind spinless atoms in a quasi ID optical lat- 
tice, the creation, a|, or annihilation operator a; in 
'creates' or 'annihilates' an atom in the Zth well of the 
optical potential, in a Wannier state ipii^x). The param- 
eter J is the hopping matrix element, the parameter W 
the on-site interaction energy for atoms sharing one and 
the same well. The Hilbert space of the system Q is 
spanned by Fock states given (in the coordinate rep- 
resentation) by symmetrised products of single-particle 
Wannier states. Assuming a finite size of the lattice, 
I = 1, . . . , L, and a finite number N of the atoms (which 
is a conserved quantity), the dimension of the Hilbert 
space isM ={N + L- 1)!/7V!(L - 1)!. 

It is worth noting that instead of the Wannier basis 
one can use a basis given by symmetrised products of 
Bloch waves 4>k{x) — {l/^/L)'^iexp{iKl)'ilji(x), where 



K = 2'!Tk/L is the single-particle quasimomentum {n = 2t: 
corresponds to the centre of the Brillouin zone). In sec- 
ond quantisation, this change of the basis amounts to the 
substitution 6/j = {l/\/L)'^^exp{i2Trkl/L)ai. Then the 
Hamiltonian Q takes the form (up to an additive term 
WN/2) 
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where d{k) = 1 if fc is a multiple of L, and d{k) — 
otherwise (conservation of total quasimomentum). 

We proceed with our analysis of the spectrum of 
the BH model. In the limits J = (no hopping) 
and W = (non- interacting bosons), the BH model is 
completely integrable, with eigenstates given by sym- 
metrised products of the single-particle Wannier states 
. . . , n/, . . . , rti), or by symmetrised products of Bloch 
states |ni, . . . , n^, . . . , n^), respectively. The correspond- 
ing eigenvalues are 



for J = 0, and 
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for W = 0. Note that in both cases the majority of 
the energy levels are multiply degenerate, i.e. different 
combinations of integers ni (uk) may result in the same 
values of {v). 

Now we demonstrate the non-integrability of the BH 
model for J,W ^ 0. The left panel in Fig. ^ shows its 
spectrum for finite size i = = 5, as a function of 
the parameter u = W , with J — 1 — u. The first step 
of the analysis is to account for the trivial symmetry of 
the spectrum, due to the translational invariance of the 
Hamiltonian. This implies that the Hamiltonian H1I2|I 
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FIG. 1: Energy level dynamics of the BH model, as a func- 
tions of the parameter u = W, J — 1 — u, for L = N = 5 and 
periodic boundary conditions. The left panel shows the com- 
plete spectrum, the right panel singles out the 'odd n — 2n^ 
symmetry component. 



factorises as 



(5) 
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where k = 27rfc/L is now the total quasimomentum. 
Thus, the spectrum in Fig. ^a) is a superposition of L 
independent spectra. Next we account for the 'odd-even' 
symmetry of the k = tt and k = 2n subspaces. The lev- 
els of the odd k — 2t: symmetry are shown in the right 
panel of Fig. ^ No crossings between levels are left (on 
the scale of the figure, some tiny avoided crossings are 
not resolved), and, thus, no further decomposition of the 
spectrum is possible. This proves the non-integrability 
of the system J^l • It is interesting to note in this context 
that the ID Hubbard model for fermions is integrable for 
arbitrary u. In the latter case, in addition to the trivial 
(it-independent) symmetries (i.e., quasimomentum and 
odd-even symmetry), the system has a number of non- 
trivial (u-dependent) integrals [1,01 • Therefore, the spec- 
trum of the ID Hubbard model can be decomposed fur- 
ther until only one (possibly degenerate) level is left. 

The above non-integrability of our system (manifest in 
the avoided crossings in the right panel of Fig. ^) does not 
yet guarantee that the spectrum of the BH model has the 
universal statistical properties associated with quantum 
chaos. Loosely speaking, the latter requires a delocali- 
sation of the eigenfunctions of the system in any generic 
basis. To quantify such delocalisation, we consider the 
Shannon entropy 



S{u) = { 



where the Cj are the expansion coefficients of an arbitrary 
eigenstate of the system in a given basis, and the angu- 
lar brackets denote an average over all eigenstates. It is 
reasonable to choose as a basis the 'I'-basis', when our 
starting point is it = (no interaction), or the '/i-basis', 
when starting at m = 1 (no hopping). First let us clarify 
the notion of v- and /i-basis. As mentioned above, and 
as apparent from Fig. ^ the spectrum of the system at 
M = 1 ( J = 0) consists of a number of degenerate lev- 
els with energies given in Eq. (O. An arbitrarily small 
J removes this degeneracy, and the levels split into a /i- 
band with a width proportional to J. Moreover, there is 
a vicinity of J = where the eigenfunctions of the sys- 
tem do not depend on J. These states form the /i-basis. 
More formally, the /i-basis is given by the eigenstates of 
the Hamiltonian ^ where the matrix elements between 
Fock states with different /i are set (by hand) to zero. 
Analogously, the i^-basis is given by the eigenstates of 
the Hamiltonian |(2J) where the matrix elements between 
states with different v are set to zero. 

The entropy (jSJ calculated in the z/-basis (solid line) 
and in the /t-basis (dashed line), respectively, is de- 
picted in the left panel of Fig. |3 for a filling factor 
n = N/L = 1. There is a relatively wide interval of u 
where the functions S'^'^^ (u) and 5*^^^ (u) simultaneously 
take large values. Thus, the eigenfunctions of the sys- 
tem are essentially delocalised in either of the two ba- 
sis sets, and one may expect universal spectral statis- 
tics. Indeed, a statistical analysis of the spectrum con- 
firms this expectation. The histogram in the upper panel 
of Fig. 121 shows the distribution of eigenenergies E for 
L = N — 8 and u = 0.3. The numerically obtained den- 
sity of states is well approximated by a Gaussian distribu- 
tion, p{E) ~ exp[— (_E — i^)^/(T^], where E and cr arc fit- 
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FIG. 2: Left panel: Mean entropy @ in the //-basis (solid 
line) and in the /i-basis (dashed line), as a function of the 
interaction parameter u {N = L = 8). Right panel: Entropy 
nun of the individual eigenstates (dots), for u = 0.3. 
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FIG. 3: Density of states (upper panel) and integrated level 
spacing distributions (lower panel) for L — N = 8 and u — 
0.3. The solid line in the upper panel is the best Gaussian 
fit to the numerically obtained density of states. The dash- 
dotted, dashed and dotted lines in the lower panel represent 
the random matrix predictions 17191 for Poissonian, Gaussian 
orthogonal (GOE), and two superimposed, independent GOE 
ensembles, respectively. The solid lines are obtained from the 
numerical data for k — 2n (which reproduce the expectation 
for superimposed GOE spectra), and for k = 2-k/8 (which 
faithfully obey the GOE prediction). 



ting parameters |0 . Using E and a from the best fit we 
unfold the spectrum and extract the nearest-neighbour 
level statistics for k — 2tt and k = 2tt/L. The result is 
presented in the lower panel of Fig. O by the integrated 
level spacing distribution I{s), s = AEp{E) (solid line), 
in comparison to the Random Matrix Theory predictions 
for the Poissonian (dash-dotted), 



/(s) = 1 - exp(-s) , 
and the Gaussian orthogonal ensemble (dashed), 



I{s) 1 - exp ( — 



(7) 
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as well as for the case resulting from a superposition of 
two independent GOE spectra (dotted line). 
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The numerical data faithfully follow the distributions 
(|8I9|) for quasimomentum k = 2n / L and k, = 2n, respec- 
tively, where the odd and the even part of the spectrum 
are superimposed in the latter case. 

We now address the conditions for quantum chaos to 
prevail in the BH model. Obviously, the above criterion 
based on the mean entropy ^ provides only a rough 
estimate on the parameter region where one may expect 



FIG. 4: Level dynamics of the modified BH model IIH for u = 
0.1, K = 2n, and L — N = 5. The grey shaded energy level 
is the diabatic approximation for the ground state evolution, 
which neglects the interaction induced avoided crossings. 



an irregular spectrum. Moreover, besides the interaction 
parameter u and the filling factor n, also the energy E 
should be considered as a relevant parameter. As a an 
illustration of this statement, the right panel in Fig. |21 
shows the entropy of the individual eigenstates. 



SiE) =inm[S^''\E),S^''\E)] , 
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which are labeled here by the individual level's energy E. 
Obviously, there is a nonvanishing fraction of localised 
eigenstates associated with the lower and upper parts of 
the spectrum. Hence, one has to distinguish regular and 
chaotic spectral components jllj . In the presently con- 
sidered case (m — 0.3, n — 1) the regular component is 
negligible (less than 1/10 of the total number of levels) 
and, hence, the level spacing distribution follows the uni- 
versal Wigner-Dyson law. However, when the ratio J/W 
approaches either one of the integrable limits, the regu- 
lar component gradually increases at the expense of the 
irregular one, what causes a deviation from the univer- 
sal distribution. The same holds true if we decrease the 
value of the filling factor. In this latter case the chaotic 
component (typically located in the central part of the 
spectrum) is small even for J/W ~ 1. Further studies 
are needed for a precise mapping of the chaotic region of 
the BH model in the parameter space spanned by u, n, 
and E. 

Let us now discuss the physical manifestations of the 
irregular spectrum of the BH model, i.e. its relevance 
for ongoing experiments with cold atoms in optical lat- 
tices la, [Tj, [l3 ■ Obviously, to probe the irregular part 
of the spectrum one has to excite the system. This can 
be done, for example, by applying a static force, i.e., 
adding a term F'Y^ilfii to the Hamiltonian Note 
that the translational symmetry of the system, broken by 
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this term, can be recovered by a suitable gauge transfor- 
mation |l5l | . This leads to a time-dependent Hamiltonian 
of the form 

Hit) (^e-^* + + Y E ^'(^'-1) ' 

(11) 

and a parametric dependence of the eigenvalues of the 
Hamiltonian on 9 = Ft, as illustrated in Fig. 01(for 
L ~ N ~ 5, K ~ 2tt, as in Fig. ^ and u = 0.1). Let 
us assume for a moment that there are no atom-atom 
interactions and, thus, all avoided crossings in Fig.^are 
true crossings. Then, in the course of time, the static 
force will drive the system along the continuation of the 
ground state (the grey shaded energy level in Fig.^J, rep- 
resenting nothing else than bona fide Bloch oscillations, 
where the system comes back to its initial state after 
one Bloch period Tg = 2-k/F. It is clear, however, that 
for It ^ the system generally does not come back to 
the initial state, because of Landau-Zener transitions at 
the avoided crossings encountered during one Bloch cycle 
[T^. In other words, after the Bloch cycle is completed, 
some excited states will be populated as a consequence 
of (a)diabatic transitions along the diabatic continuation 



of the ground state. (The branching ratio between di- 
abatic and adiabatic transitions is obviously determined 
by the time derivative of 6, i.e., by the field amplitude F.) 
In the case of strong static forcing, F J, this excita- 
tion process is reversible, leading to quasiperiodic Bloch 
oscillations 17] . In contrast, for moderate forcing (and 
chaotic instantaneous spectrum T^) the excitation is ir- 
reversible, leading to decoherence of the one-particle den- 
sity matrix and, as a consequence, to decay of the Bloch 
oscillations jl3| • Since Bloch oscillations of the atoms are 
easily measured in laboratory experiments |l2l lll^ , their 
irreversible decay provides a direct indication of chaos in 
the BH system [13^ . 

In conclusion, we have shown that the ID Bose- 
Hubbard model is a generic quantum chaotic system with 
characteristic spectral irregularity. This is in stark con- 
trast to the ID Hubbard model for fermions, which is 
known to be an integrable system with regular spectrum. 
It is also worth noting that the many-site character of the 
BH model is crucial for the appearance of quantum chaos 
- the two site BH model, which is often used to study 
Bose Einstein condensates of cold atoms in double-well 
potentials Fh? , has a regular spectrum. 
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